clear; close all; clc;

rng(42);

n_samples = 100;

wavelength = 400 + 200 * rand(n_samples, 1);
refractive_index = 1.45 + 0.3 * rand(n_samples, 1);
transmittance = 0.85 + 0.1 * rand(n_samples, 1);
thickness = 1 + 4 * rand(n_samples, 1);
curvature_radius = 50 + 100 * rand(n_samples, 1);

transmittance = transmittance - 0.05 * (thickness - 3);

X = [wavelength, refractive_index, transmittance, thickness, curvature_radius];
feature_names = {'波长(nm)', '折射率', '透过率', '厚度(mm)', '曲率半径(mm)'};

X_standardized = zscore(X);

[coeff, score, latent, tsquared, explained] = pca(X_standardized);

fprintf('=== 光学元件参数PCA分析结果 ===\n\n');
fprintf('数据集大小: %d个样本 × %d个特征\n\n', n_samples, size(X,2));

fprintf('主成分贡献率:\n');
for i = 1:length(explained)
    fprintf('PC%d: %.2f%%\n', i, explained(i));
end
fprintf('前两个主成分累计贡献率: %.2f%%\n\n', sum(explained(1:2)));

fprintf('主成分载荷矩阵:\n');
fprintf('特征\t\tPC1\t\tPC2\t\tPC3\t\tPC4\t\tPC5\n');
for i = 1:length(feature_names)
    fprintf('%s\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\n', ...
        feature_names{i}, coeff(i,1), coeff(i,2), coeff(i,3), coeff(i,4), coeff(i,5));
end

figure('Position', [100, 100, 1200, 800]);

subplot(2,3,1);
pareto(explained);
xlabel('主成分');
ylabel('贡献率 (%)');
title('主成分贡献率');
grid on;

subplot(2,3,2);
scatter(score(:,1), score(:,2), 50, 'filled', 'MarkerFaceAlpha', 0.6);
xlabel(['PC1 (', num2str(explained(1), '%.1f'), '%)']);
ylabel(['PC2 (', num2str(explained(2), '%.1f'), '%)']);
title('PCA二维投影');
grid on;

subplot(2,3,3);
scatter3(score(:,1), score(:,2), score(:,3), 50, 'filled', 'MarkerFaceAlpha', 0.6);
xlabel(['PC1 (', num2str(explained(1), '%.1f'), '%)']);
ylabel(['PC2 (', num2str(explained(2), '%.1f'), '%)']);
zlabel(['PC3 (', num2str(explained(3), '%.1f'), '%)']);
title('PCA三维投影');
grid on;

subplot(2,3,4);
imagesc(abs(coeff));
colorbar;
xticks(1:5);
xticklabels({'PC1', 'PC2', 'PC3', 'PC4', 'PC5'});
yticks(1:5);
yticklabels(feature_names);
title('主成分载荷绝对值热图');
xlabel('主成分');
ylabel('原始特征');

subplot(2,3,5);
biplot(coeff(:,1:2), 'Scores', score(:,1:2), 'VarLabels', feature_names);
xlabel(['PC1 (', num2str(explained(1), '%.1f'), '%)']);
ylabel(['PC2 (', num2str(explained(2), '%.1f'), '%)']);
title('PCA双标图');

subplot(2,3,6);
corr_matrix = corr(X);
imagesc(corr_matrix);
colorbar;
caxis([-1, 1]);
xticks(1:5);
yticks(1:5);
xticklabels(feature_names);
yticklabels(feature_names);
title('原始特征相关性矩阵');

for i = 1:size(corr_matrix,1)
    for j = 1:size(corr_matrix,2)
        text(j, i, num2str(corr_matrix(i,j), '%.2f'), ...
            'HorizontalAlignment', 'center', 'FontSize', 8);
    end
end

X_reduced = score(:,1:2);

fprintf('\n=== 降维结果分析 ===\n');
fprintf('原始数据维度: %d\n', size(X,2));
fprintf('降维后数据维度: %d\n', size(X_reduced,2));
fprintf('信息保留率: %.2f%%\n', sum(explained(1:2)));

X_reconstructed = score(:,1:2) * coeff(:,1:2)' .* std(X) + mean(X);
reconstruction_error = mean(mean((X - X_reconstructed).^2));
fprintf('重构均方误差: %.4f\n', reconstruction_error);

figure('Position', [100, 100, 800, 600]);

k = 3;
[idx, centroids] = kmeans(X_reduced, k);

gscatter(score(:,1), score(:,2), idx, 'rgb', 'o', 8);
hold on;
plot(centroids(:,1), centroids(:,2), 'kx', 'MarkerSize', 12, 'LineWidth', 3);
xlabel(['PC1 (', num2str(explained(1), '%.1f'), '%)']);
ylabel(['PC2 (', num2str(explained(2), '%.1f'), '%)']);
title('基于PCA的样本聚类');
legend('类别1', '类别2', '类别3', '聚类中心');
grid on;

results_table = array2table([X, X_reduced], ...
    'VariableNames', [feature_names, {'PC1', 'PC2'}]);
writetable(results_table, 'optical_pca_results.csv');

fprintf('\n分析完成！结果已保存到 optical_pca_results.csv\n');
